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ABSTRACT 


A  computer  algorithm  was  developed  to  determine  if  an  acoustic  transmitter  can 
be  localized  based  on  estimates  of  local  angles  of  arrival  of  acoustic  signals  incident 
upon  a  receive  planar  sonar  array,  knowledge  of  the  deterministic  effects  of  the  ocean 
on  sound  propagation,  and  local  sound-speed  profiles  of  the  ocean.  The  algorithm  was 
designed  to  determine  azimuthal  and  elevation/depression  angles  to  the  transmitter  as 
well  as  computing  the  depth,  range,  cross  range,  and  line-of-sight  range  separations 
between  the  transmitter  and  the  receive  array.  The  algorithm  utilizes  ray  acoustics  and 
model-based  phase  weights  to  determine  the  transmitter's  location  relative  to  the 
receive  array's  position.  As  written,  the  algorithm  is  capable  of  solving  localization 
problems  in  which  the  transmitter  and  receiver  are  in  the  same  gradient  of  the  local 
sound-speed  profile,  provided  that  the  range  from  transmitter  to  receiver  is  not  so  great 
that  the  acoustic  signal  passes  through  a  turning  point  prior  to  reaching  the  receive 
array.  The  results  indicate  that  the  method  proposed  is  viable  for  the  class  of  problems 
for  which  it  was  designed,  and  accuracies  on  the  order  of  0.1  meters  are  obtained  for 
line-of-sight  ranges  on  the  order  of  several  kilometers.  The  angles  calculated  by  the 
algorithm  are  all  accurate  to  within  0.005  degrees. 
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I.  INTRODUCTION 

This  thesis  constitutes  one  part  of  a  long  range  project  to  develop  new  sonar 
signal  processing  algorithms  capable  of  rapidly  solving  sonar  localization  problems.  At 
present,  the  solution  of  the  sonar  fire  control  problem  can  require  a  considerably 
longer  time  than  that  required  for  most  other  types  of  fire  control  problems.  The  long 
time  required  to  achieve  a  solution  can  cause  a  significant  degradation  in  a  ship's 
ability  to  avoid  counterdetection,  due  to  continuously  decreasing  range  to  the  target 
during  problem  solution.  A  sonar  system  capable  of  rapid  target  localization  without 
requiring  own  ship's  maneuvers  would  greatly  enhance  the  capabilities  of  our  ships, 
and  allow  for  weapon  firings  at  longer  ranges. 

The  research  question  investigated  in  this  thesis  is  whether  or  not  it  is  possible 
to  develop  an  algorithm  which  utilizes  estimates  of  the  local  angles  of  arrival  of 
acoustic  signals  incident  upon  a  planar  sonar  array,  knowledge  of  the  deterministic 
effects  of  the  ocean  medium  on  sound  propagation,  and  local  sound-speed  profiles  of 
the  ocean,  to  locate  an  acoustic  transmitter,  both  in  azimuthal  angle  and 
elevation/depression  angle.  In  addition  the  model-based  localization  algorithm 
(hereafter  referred  to  as  the  'localization  algorithm')  was  designed  to  provide  the  range, 
depth,  cross  range,  and  line-of-sight  range  between  the  acoustic  transmitter  and  the 
receive  array. 

Ray  acoustics  provides  methods  of  determining  ranges  and  propagation  angles 
for  transmission  of  acoustic  signals  in  inhomogeneous  media  [Ref.  1:  sect.  6.2].  The 
deterministic  effects  of  the  inhomogeneous  ocean  medium  on  acoustic  signals  are  well 
known.  From  a  transmitter  in  a  known  position,  it  is  possible  to  develop  ray  traces 
that  illustrate  the  propagation  of  acoustic  signals  through  the  ocean  medium.  The 
intent  here  is  to  use  this  knowledge  of  sound  propagation  to  find  the  transmitter's 
location  based  on  the  estimated  angles  of  arrival  at  a  receive  array.  The  estimates  of 
the  local  angles  of  arrival  are  obtained  from  a  frequency  domain  adaptive  beamforming 
algorithm  developed  by  Ziomek  and  Chan  [Ref.  2].  This  algorithm  performs  frequency 
domain  adaptive  beamforming  for  planar  sonar  arrays  using  a  modified  complex  LMS 
adaptive  algorithm.  The  algorithm  generates  estimates  of  the  local  angles  of  arrival, 
namely,  the  azimuthal  and  elevation  depression  angles,  of  incoming  acoustic  signals. 


However,  in  a  real  ocean  environment,  these  local  angles  of  arrival  do  not  reflect  the 
true  line-of-sight  angles  to  the  target. 

The  localization  algorithm  uses  the  angle-of-arrival  estimates,  plus  typical 
sound-speed  profiles  that  are  normally  available  to  ships.  In  addition,  it  was  found 
that  one  more  piece  of  information  is  required  to  localize  the  target.  This  information 
is  a  model-based  phase  weight  which  is  part  of  a  model-based  signal  processing 
algorithm  developed  by  Ziomek  and  Blount  [Ref.  3].  These  phase  weights  are  used  to 
"correct  for  deterministic,  ocean  medium,  phase  effects  due  to  ray  bending  as  a  signal 
propagates  in  the  inhomogeneous  ocean  medium  whose  index  of  refraction  (sound- 
speed  profile)  is  a  function  of  depth."  [Ref.  3]  The  phase  weights  were  originally 
developed  as  part  of  an  underwater  acoustic  communication  problem  in  which  receiver 
and  transmitter  locations  were  known.  The  form  of  the  phase  weights  used  will  be 
presented  in  Chapter  II. 

For  the  problem  investigated  in  this  thesis,  transmitter  location  is  unknown  a 
priori  and,  therefore,  the  model-based  phase  weights  cannot  be  determined  in  exactly 
the  same  manner  as  was  done  in  the  algorithm  developed  by  Ziomek  and  Blount 
[Ref.  3].  The  usefulness  of  the  localization  algorithm  developed  in  this  thesis  is  based 
on  the  availability  of  the  model-based  phase  weights.  The  research  done  here  is  a 
feasibility  study  of  the  ability  to  localize  an  acoustic  transmitter  if  the  phase  weights 
were  available.  The  development  of  an  algorithm  to  generate  the  model-based  phase 
weights  was  not  the  subject  of  this  research. 

The  localization  algorithm  is  limited  to  solving  a  particular  class  of  problems. 
The  localization  algorithm  is  designed  to  accommodate  vertical  variations  in  sound- 
speed  profile  or,  sound-speed  profiles  that  are  functions  of  depth  only.  Horizontal  or 
range  variations  in  sound-speed  profile  were  not  examined  in  this  initial  study  because 
they  constitute  only  a  relatively  small  portion  of  ocean  areas.  Additionally,  the 
transmitter  and  receiver  are  assumed  to  both  be  within  the  same  sound-speed  gradient. 
Finally,  all  case  studies  were  conducted  based  on  the  assumption  that  the  receiver  was 
in  close  enough  proximity  to  the  transmitter  so  that  the  acoustic  signal  had  not  passed 
through  a  turning  point  prior  to  reaching  the  receiver.  A  turning  point  is  defined  as 
the  point  along  a  ray  path  at  which  the  angle  of  propagation  is  90  degrees  with  respect 
to  the  positive  Y,  or  depth,  axis.  These  three  restrictions  were  necessary  to  limit  the 
scope  of  the  initial  study  to  a  size  that  would  allow  for  a  complete  verification  of  the 
localization  technique  proposed,  in  the  time  allotted  for  the  study. 
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Chapter  II  describes  the  theory  used  to  develop  the  localization  algorithm.  An 
overview  of  the  problem  and  its  geometry  is  presented,  and  then  the  computations 
leading  to  the  algorithm  are  discussed.  Finally,  the  limitations  of  the  algorithm  are 
presented. 

Chapter  III  consists  of  the  computer  simulation  results  and  an  explanation  of  the 
implementation  of  the  theory  in  a  computer  algorithm.  The  output  from  the 
localization  algorithm  is  compared  to  the  known  geometry,  and  a  comparison  of 
double  precision  versus  single  precision  results  is  included.  Additionally,  the  program 
is  investigated  to  determine  if  errors  develop  as  a  function  of  the  transmission  angle 
and  or  depth  separation.  As  will  be  shown  in  Chapter  II,  the  roots  of  a  fourth-order 
polynomial  must  be  determined  to  find  the  angle  of  transmission  at  the  source.  The 
roots  for  the  fourth-order  polynomial  are  found  through  use  of  an  International 
Mathematical  Subroutine  Library  (IMSL)  subroutine  and  are  verified  by  comparison 
with  graphs  of  the  function.  These  graphs  also  assist  in  determining  the  correct  root  to 
use  during  problem  solution. 

In  Chapter  IV,  conclusions  and  recommendations  are  presented. 
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II.  THEORY 

A.       PROBLEM  OVERVIEW  AND  GEOMETRY 

Traditionally,  the  localization  of  acoustic  transmitters  by  ships  has  been  carried 
out  by  obtaining  many  lines  of  bearing  to  the  transmitter,  and  comparing  these  with 
own  ship's  motion  to  develop  a  geographic  picture  of  the  transmitter's  motion.  This 
method  is  time  consuming  and  usually  very  lacking  in  terms  of  accuracy.  Due  to  the 
nature  of  the  deterministic  effects  of  the  ocean  medium,  a  great  deal  of  information  is 
contained  in  the  angles  at  which  acoustic  energy  arrives  at  the  receiver.  Extraction  of 
this  information  from  the  local  angles  of  arrival,  while  not  a  simple  task  in  of  itself, 
would  greatly  simplify  the  problem  of  target  localization. 

As  a  first  step  in  exploiting  the  information  contained  within  the  local  angles  of 
arrival,  a  geometry  must  be  assumed  for  the  problem.  Figure  2.1  illustrates  the  general 
three-dimensional  geometry  used  in  the  development  of  the  method  of  target 
localization  presented  here. 

From  Figure  2.1  the  following  definitions  are  apparent: 

•  Xq,  yQ,  z0  rectangular  coordinates  of  the  transmitter  in  meters. 

•  xR,  yR,  zR  rectangular  coordinates  of  the  center  of  the  receive  planar 

array  in  meters. 

•  AX,  AY,  AZ  cross  range,  depth,  and  Z  coordinate  separations,  respectively, 

in   meters   between   the   transmitter   and   the   receive   array, 
where: 

AX  =  xR  -  x0 

AY  =  yR  -  y0 

AZ  =  zR  -  z0 

•  AR  polar  radial  distance  in  meters  from  the  transmitter  to  the 

receive  array. 

Note:  AR2  =  AX2  +  AZ2 

•  HDLTR  polar  radial  distance  in  meters  that  a  ray  would  travel  in  a 

homogeneous  medium  (constant  sound-speed  profile)  between 
depths  yQ  and  yR  based  on  an  angle  of  transmission  of  p(y0). 

•  HDLTX,  HDLTZ    distances  in  the  X  and  Z  directions,  respectively,  that  a  ray 

would  travel  in  a  homogeneous  medium  between  depths  yQ 
and  yn  based  on  an  angle  of  transmission  of  p(yQ). 

•  FIRLOS  line-of-sight  range  that  a  ray  would  travel  in  a  homogeneous 

medium  between  depths  y0  and  yR  based  on  an  angle   of 
transmission  of  P(y„). 

12 
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Figure  2.1     System  Geometry. 
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Note:   HRLOS2  =  HDLTX2  +  AY2  +  HDLTZ2. 

•  RLOS  line-of-sight  range  between  the  transmitter  and  the  center  of 

the  receive  array. 

Note:    RLOS2  =  AX2  +  AY2  +  AZ2 

•  P(y0)  initial  angle  of  propagation  (angle  of  transmission),  measured 

with  respect  to  the  positive  Y  axis,  of  the  acoustic  signal  at 
source  depth  yQ  meters. 

•  P(Yr)  angle   of  arrival   of  incident  plane  wave   field  at  depth  yrj^ 

meters. 

•  PLOS  the  line-of-sight  angle,  as  measured  from  the  positive  Y  axis, 

between  the  transmitter  and  the  receive  array. 

Note  that  in  Figure  2.1  the  positive  Y  axis  is  defined  in  the  direction  of  increasing 
depth,  or  in  the  downward  direction.  The  coordinate  system  shown  in  Figure  2.1  is 
applicable  for  any  relative  positioning  of  the  transmitter  and  receive  array,  even  if  AX, 
AY,  and/or  AZ  are  negative.  Thus,  the  algorithm  will  work  for  any  direction  of  arrival 
of  the  incident  acoustic  plane-wave  field. 

The  receive  array  is  assumed  to  possess  knowledge  of  its  own  depth.  In  addition, 
the  receive  array  will  have  available  estimates  of  arrival  direction  cosines  associated 
with  the  local  angles  of  arrival.  These  estimates  are  computed  by  the  frequency 
domain  adaptive  beamforming  algorithm.  From  these  known  quantities  and 
information  about  the  local  sound-speed  profile,  the  transmitter's  location  with  respect 
to  the  receiver  shall  be  determined. 

B.   TRANSMITTER  LOCALIZATION  THEORY 

Energy,  whether  it  is  acoustic  or  electromagnetic,  will  refract  as  it  passes  from  a 
medium  with  index  of  refraction  nt  into  a  medium  with  index  of  refraction  n2,  provided 
that  n}  *  nr  In  this  study,  the  ocean  volume  is  characterized  by  a  one-dimensional 
index  of  refraction  (sound-speed  profile)  that  is  a  function  of  depth.  Snell's  law  is 
given  by  [Ref  1:  p.  218], 

sin  P(y)  _  sin  P(yQ) 
c(y)  c(yQ) 

where  c(y)  is  the  speed  of  sound  in  meters  per  second  at  a  depth  y.    From  Snell's  law  a 
ray  parameter  may  be  defined  as 
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sin  P(y0)  __  sin  p(yR)  =  sin  P(yTP)  1 

c(y0)  c(yR)  c(yTp)         c(yTp) 

where: 

•  b  is  the  ray  parameter. 

•  yjp  is  the  depth  of  a  turning  point.    A  turning  point  is  defined  as  the  point 

along  a  ray  path  at  which  the  angle  of  propagation,  P(y-j-p),  is  equal  to 
90  degrees. 

At  this  point  P(yR)  is  known,  since  the  direction  cosine 

v(yR)  =  cos  P(yR)  (2.3) 

is  calculated  by  the  frequency  domain  adaptive  beamforming  algorithm.  The  speed  of 
sound  at  depth  yR,  denoted  c(yR),  is  normally  known  aboard  ship  as  a  result  of 
measurements  made  by  onboard  sonar  systems. 

It  is  assumed  that  the  sound-speed  profile  is  a  linear  function  of  depth  with 
constant  gradient.  In  most  areas  of  the  ocean  this  is  a  good  approximation  if  both  the 
transmitter  and  the  receive  array  are  in  the  same  portion  of  the  sound-speed  profile.  A 
typical  sound-speed  profile  is  shown  in  Figure  2.2.  The  parameter  g  is  the  constant 
gradient  of  the  sound-speed  profile  in  seconds"1.  From  the  surface  to  about  100  meters 
a  positive  gradient  is  typically  observed  with  a  gradient  g  ~  +0.016  sec 
[Ref.  4:  p.  30],  [Ref.  5:  p.  401].  Below  100  meters  a  negative  gradient  is  present,  and  in 
this  example  g  «  -0.02956  sec"1.  Finally,  at  depths  between  700  to  1500  meters 
[Ref.  4:  p.  32]  the  gradient  reverts  to  a  positive  value  of  g  *  +0.017  sec"1 
[Ref.  5:  p.  401].  The  value  of  g  in  the  negative  portion  of  the  gradient  was  computed 
by  assuming  the  speed  of  sound  to  be  1500  meters  per  second  at  the  ocean  surface  and 
1475  meters  per  second  at  a  depth  of  1000  meters  [Ref.  6:  p.  3].  A  depth  of  1000 
meters  was  chosen  as  the  starting  point  of  the  second  positive  gradient.  The  negative 
gradient  was  then  calculated  to  fit  between  the  positive  gradients.  Based  on  the 
assumption  that  both  the  transmitter  and  the  receive  array  are  in  the  same  gradient  of 
the  sound-speed  profile,  the  speed  of  sound  at  depth  y  can  be  found  from 

c(y)  =  c(yn)  +  g(y  -  yn)  .  (2.4) 
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Speed  of  Sound 


Q. 

CD 

Q 


Range 


g  =  0.016/sec 

N.  Sound-Speed    Profile 
\   100  m 

g  =  -0.0296/sec 

- ► 

1000  m  ^ 

\    g  =  0.017/sec 

Figure  2.2    Typical  Sound-Speed  Profile. 

The  radius  of  curvature  that  describes  the  arc  of  the  circle  followed  by  an 
acoustic  field  propagating  through  this  medium  is  then  [Rcf.  1:  p.  237] 


R    =        C(>Q)         =        C(VR} 
c       |g  sin  p(yn)|        |gsinP(yR)| 


(2.5) 


All  the  terms  on  the  far  righthand  side  of  Equation  2.5  are  known.    Equation  2.4  may 
be  rewritten  as  expressed  by  Ziomck  [Rcf.  1:  p.  238] 


V  "  >'n  + 


c(y)  -  c(y0) 


(2.6) 


Therefore, 
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1  1 

AY  =  yR  -  yQ  =  —  c(yR)  -  —  c(y0)  (2.7) 


and,  from  equations  2.1  and  2.2, 
sin  P(yR) 


c 


(yR)  = -*-  (2.8) 


and 


sin  p(yn) 
c(y0)  -  — ^  •  (2.9) 


Combining  equations  2.8  and  2.9  with  equation  2.7  it  is  readily  observed  that, 


AY  =  >'R  -  >'o  =  7-  sin  P(>'R>  -  7-  sin  P(y0>  <2-10) 

bg  bg 


or, 

AY  =  >'R  ■  >'o  =  a  sin  P(>R)  _  a  sin  P<V  (2-u) 

where 


be 


The  only  unknowns  now  in  equation  2.1 1  are  P(y0)  and  AY.   Also  note  that 

R    =  |a|  =  radius  of  curvature.  (2.13) 

The  radial  distance  AR  shown  in  Figure  2.1   can  be  found  by  utilizing  the 
following  equation  [Ref.  1:  p.  238]: 
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z  =  zo+    ^IV  ; [cos  P(V  - cos  Wy)]  (2-14) 

g  sinP  (yQ) 

which  is  the  Z  coordinate  of  a  ray  propagating  in  the  YZ  plane.  In  this  thesis  a  more 
general  class  of  problem  is  assumed  so  that  the  coordinate  axes  can  remain  fixed 
relative  to  the  platform  on  which  the  planar  array  is  mounted.  Therefore,  in  a  three- 
dimensional  system,  z  and  zQ  are  replaced  by  the  polar  coordinates  r  and  rQ  to  give 

r  =  ro  +     c(o"}  ;  tcos  P(y0)  - cos  P(y)i .  (2-15) 

g  sinp  (yn) 


and,  as  a  result, 


AR  =  rR  -  rQ  =  —  [cos  p(yQ)  -  cos  0(yR)]  (2.16) 


or 


AR  =  rR  -  r0  =  a  cos  p(yQ)  -  a  cos  p(yR)  .  (2.17) 

The  only  unknowns  in  equation  2.17  are  P(y0)  and  AR.  Also,  note  in  Figure  2.1  that  if 

AX  =  0,  (2.18) 

then 

AR  =  AZ.  (2.19) 

At  this  point  ray  acoustics  cannot  provide  any  further  information  to  develop  a 
solution  to  the  problem.  However,  a  model-based  phase  weight  for  a  planar  sonar 
array,  similar  to  that  shown  in  Figure  2.3  ,  can  be  used  to  localize  the  transmitter.  As 
derived  by  Ziomek  and  Blount  [Ref.  3] 

4>n(D  =  -27tfYndY  +  OM0(f,n)  n  =  -(N-l)/2,...,0,...,(N-l)/2  (2.20) 
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m  =  2 


Figure  2.3     Receive  Planar  Arrav  Gcomctrv. 


19 


where 

-vRf 
fy— A-  (2.21) 

c(yR) 

vB  =  cos  P(yT)  ,  (2.22) 

and: 

•  <D  (f)  is  the  phase  weight  in  the  Y  direction  associated  with  element  (m,n) 

in  the  receive  array. 

•  f  is  the  frequency  of  the  transmitted  electrical  signal. 

•  <I>iyjr)(f,n)    is  the  model-based  phase  weight  which  is  related  to  the  deterministic 

angle  modulation  performed  by  the  ocean  medium  on  the  transmitted 
electrical  signal  as  a  function  of  depth  [Ref.  3]. 

•  yj  is  the  depth  of  the  transmit  array. 

•  dy  is  the  interelement  spacing  in  the  Y  direction  associated  with  the 

receive  array. 

Equation  2.20  describes  the  phase  weights  in  the  Y  direction  that  a  planar  sonar 
array  using  the  three-dimensional  FFT  beamformer  presented  by  Ziomek  and  Blount 
[Ref.  3]  would  use  to  receive  an  acoustic  signal  transmitted  from  a  depth  yj  and 
received  at  a  depth  yn.  This  equation  can  be  seen  to  consist  of  two  parts.  The  first 
portion  is  the  term  ^TtfyndY  which  is  the  phase  weight  used  in  traditional  beam 
steering.  The  second  part,  <Pjyrr\(f,n),  is  further  described  by  Ziomek  and  Blount 
[Ref.  3]  as 

<I>MD(f,n)  =  -  [k(yT)/2vB]{[c(yT)/g][n'D(yR  +  ndy)  -  1]  +  AYn}  (2.23) 

where  the  wave  number  in  radians  per  meter  as  a  function  of  depth  yj  is 

k(yT)  =  2Jlf/c(yT)  ,  (2.24) 


f  =  f '  +  kfn  ,  k  =  -K,. ..,(),..-., K  ,  (2.25) 


'o 
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n  D(yR  +  ndY)  =         '  ?T  ,  (2.26) 

c(yT)  +  g\\n 

and, 

AYn  =  >'R  -  >'T  +  ndY  •  <2-27) 

where: 

•  f     is  the  carrier  frequency  in  hertz, 

•  fQ     is  the  fundamental  frequency  in  hertz  of  the  finite  Fourier  series  representation 

of  the  complex  envelope  of  the  transmitted  electrical  signal,  and 

•  K     is  the  highest  harmonic  used  in  the  finite  Fourier  series. 

The  term  n'^  defines  an  index  of  refraction  which  is  corrected  for  the  distance 
that  the  (m,n)  element  in  the  receive  array  is  offset  in  the  Y  direction  from  the  center 
of  the  array.  This  compensation  is  provided  by  AY  ,  which  computes  the  depth 
separation  between  the  center  of  the  transmit  array  and  the  element  (m,n). 

When  using  these  model-based  phase  weights  it  should  be  noted  that  yj  is 
equivalent  to  yQ  of  Figure  2.1.   Additionally, 

vB  =  cos  P(yT)  -  cos  P(y0)  .  (2.28) 

Dividing  equation  2.24  by  2vg  yields 

k(vT)  27tf         1  71  f 

-^-  = =  .  (2.29) 

2vB         c(yT)    2vB       c(yT)vB 

The  term  n'n(yn  +  ndy)  -  1  may  also  be  rewritten  as 

nD(vR  +  ndY)-,%(^AY/,  (2.30, 

and,  as  a  result, 
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c(vT)  -  c(vT)  -  gAY„ 
n  D(yR  +  ndy)  -  1  =       '?         'J      *      "  (2.31) 

u  c(yT)  +  gAYn 

n  D(yR  +  ndY)  -  1  =  ~*J"  •  (2-32) 

Therefore,  substituting  equations  2.32  and  2.29  into  equation  2.23  yields 

-Ttf        [-c(vT)AY    +  c(yT)AY    +  gAY  2] 
0>Mn(f,n)  =  {—^- 2 ^- 2 2_L  (2.33) 

c(>T>vB  W>T>  +  §AYJ 

-Ttf  eAY  2 

c(yT)vB   [c(yT)  +  gAYJ 

Expanding  the  denominator  of  the  second  term  on  the  right  side  of  equation  2.34 
results  in 

c(yT)  +  gAYn  =  c(yT)  +  g(yR  -  yT  +  ndy)  •  (2.35) 


From  equation  2.4  it  can  be  seen  that 

c(yy)  +  g(yR  -  yj  +  ndy)  =  c(yR  +  ndy)  •  (2.36) 

Therefore, 

c(yT)  +  gAYn  =  c(yR  +  ndy)  .  (2.37) 

Substituting  equation  2.37  and  equation  2.22  into  equation  2.34  gives 

-TtfeAY  2 

<Dvm(f,n)  = 2 ,  (2.38) 

MD  c(yT)c(yR  +  ndy)  cos  P(yT) 

From  equation  2.2,  with  y-y  —  y0 
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sin  B(vT) 

c(yT)  = r-^-  (2.39) 

b 


and,  as  a  result, 


-TTfbeAY  2 

Ovm(r,n)  =  — -2 .  (2.40) 

c(vR  +ndy)  sin  p(vT)  cos  p(yT) 


From  equation  2.12, 


a  =  —  (2.12) 

bs 


or 


1 

—  =  bg  .  (2.41) 

a 


Using  equation  2.41  in  equation  2.40  yields 


-TlfAY   2 

Ovin(f,n)  =  J (2.42) 

ac(yR  +  ndy)  sin  p(y-j-)  cos  P(yy) 


where 


AYn  =  yR  -  yT  +  ndy  =  AY  +  ndy  •  (2.43) 

If  the  center  element  of  the  receive  array  is  chosen  as  the  element  at  which  the 
phase  weight  0>^j-)(f,n)  is  calculated,  then  n  =  0,  and 

AYQ  =  AY  =  yR  -  yT  .  (2.44) 

Therefore,  at  n  =  0 
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-TtfAY2 
mu  ac(yR)  sin  P(yT)  cos  p(yT) 

Now  squaring  both  sides  of  equation  2.10  will  result  in 

AY2  =  a2sin2p(y0)  -  2a2sin  p(yR)  sin  p(y0)  +  a2sin2p(yR)  (2.46) 

which  can  be  used  in  equation  2.45  to  replace  AY2.   Now  let 

x  =  sin  P(y0)  (2.47) 

and 

y  =  cos  P(y0)  .  (2.48) 

The  x  and  y  defined  in  equations  2.47  and  2.48  are  not  the  x  and  y  coordinates 
related  to  Figure  2.1.  Rather,  this  x  and  y  are  merely  dummy  variables  to  be  used  in 
the  solution  of  equation  2.45.  Due  to  the  definitions  of  equations  2.47  and  2.48  a 
relationship  between  the  variables  x  and  y  is  apparent,  that  is, 

x2  +  y2  =  1  (2.49) 

and,  therefore, 

y  =   ±  (1  -x2)1/2.  (2.50) 

Next,  replace  yj  with  yQ  in  equation  2.45,  substitute  equation  2.46  into  equation 
2.45,  and  multiply  equation  2.45  by  ac(yR)xy.   This  results  in 

ac(yR)OMD(f,0)xy  -  -itf  [a2x2  -  2a2sin  P(yR)x  +  a2sin2P(yR)]  .  (2.51) 

Divide  both  sides  of  equation  2.51  by  Ttfa2  to  get 
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^jjk-  ^MD(f'0)x>'  =  -x2  +  2  sin  P(yR)x  -  sm2P(yR)  .  (2.52) 


Rewriting  equation  2.52  yields 


x2  +  -^■)OMD(f,0)xy  -  2  sin  P(yR)x  +  sin2p(yR)  =  0  (2.53) 


or 


Ax2  +  Bxy  +  Cx  +  D  =  0  (2.54) 


where 


A  =   1.0  ,  (2.55) 


B  -  ^^MD(f,0),  (2.56) 


C  =  -2  sin  P(yR)  ,  (2.57) 

and 

D  =  sin2  P(yR)  .  (2.58) 

Substituting  equation  2.50  into  equation  2.54  yields 

Ax2  ±  Bx(l  -  x2)1'2  +  Cx  +  D  =  0  (2.59) 

or, 

Ax2  +  Cx  +  D  =   ±  Bx(l  -  x2)1'2  .  (2.60) 
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Squaring  both  sides  of  equation  2.60  gives 

(Ax2  +  Cx  +D)2  =  B2x2(l  -  x2)  =  B2x2  -  B4x4  .  (2.61) 

Expanding  equation  2.61  yields 

A2x4  +  2ACx3  +  (2AD  +  C2)x2  +  2CDx  +  D2  =  B2x2  -  B2x4  (2.62) 


or 


(A2  +  B2)x4  +  2ACx3  +  (2AD  -   B2  +C2)x2  +  2CDx  +  D2  =  0  .  (2.63) 

To  find  the  unknown,  x,  the  roots  of  equation  2.63  must  be  computed.  These 
roots  will  also  be  the  roots  of  equation  2.59.  In  the  computer  algorithm  written  to 
implement  this  theory,  the  value 

F(x,y)  =  Ax2  ±  Bx(l  -  x2)1/2  +  Cx  +  D  (2.64) 

was  also  calculated  to  verify  the  validity  of  the  roots  found  for  equation  2.63. 
Recalling  that  equation  2.47  defined 

x  =  sin  P(yQ)  (2.47) 

and  equation  2.48  defined 

y  =  cos  P(yQ)  ,  (2.48) 

we  see  that  once  x  and  y  are  known  they  may  be  substituted  into  equations  2. 1 1  and 
2.17  to  solve  for  AY  and  AR  (since  P(}r),  the  radius  of  curvature  (a),  the  receive  array 
depth,  and  the  local  sound-speed  profile  are  all  known). 

At  this  point  AY,  AR  and  P(y0)  are  known.  The  next  values  to  be  found  are  AX, 
AZ,  RLOS,  and  PLOS.  Using  the  definitions  of  the  direction  cosines  as  presented  by 
Ziomek[Ref.  1:  p.226] 

v  =  cos  p(y),  (2.65) 
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u  =  cos  a(y),  (2.66) 

and 

w  =  cos  y(y)  (2.67) 

where: 

•  a(y)  is  the  angle  at  a  depth  y  measured  from  the  positive  X  axis  to  the 

vector  of  interest. 

•  y(y)  is  the  angle  at  a  depth  y  measured  from  the  positive  Z  axis  to  the 

vector  of  interest. 

Referring  to   Figure  2.1,  the  direction  cosine  v(y)  at  the  transmitter  depth  can  be 
written  as 


AY 

v(yn)  =  cos  P(yn)  =  (2.68) 

0  °         HRLOS 


and,  as  a  result, 

AY 

HRLOS  =  .  (2.69) 

Also  from  Figure  2.1  it  can  be  observed  that 

HRLOS2  =  HDLTR2  +  AY2  .  (2.70) 

In  ray  acoustics,  as  presented  by  Ziomek  [Ref.  1:  p. 223],  the  propagation  vector 
is  defined  as 

k  =  kxx  +  kyy  +kzz  (2.71) 

where 

kx  =  kQu  ,  (2.72) 
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kY  =  k0v  ,  (2.73) 


kz  =  k0w  ,  (2.74) 

and 

2;rf 

k0  -  — -  •  (2-75) 

c(y0) 

Therefore,  at  the  transmitter, 


kXT  =  ^7  U<>'T>  •  <2'76> 

and,  at  the  receive  array, 


kxR  =  -Jt-  u(yR)  •  (2-77) 

Additionally,  for  an  inhomogeneous  medium  which  has  a  sound-speed  profile  that  is  a 
function  of  depth  only,  it  is  known  that  [Ref.  1:  p. 223] 

kXR  =  kXT  =  constant  •  (2.78) 

Therefore,  from  equations  2.76  and  2.77, 

2?tf  2rcf 

— —  u(yR)  =  — —  u(yT)  (2.79) 

c(yR)      ^        c(yT) 

so  that 


c(Vn) 

u(>'R>  =  T^T  U(>'T)  '  (2-80> 

c(>t) 
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or, 


u(>t)  =  *n u(yR) '  (2'81) 

In  equation  2.81  u(yR)  is  supplied  by  the  beamformer,  c(yn)  is  known  by  own 
ship,  and  since  equation  2.11  has  been  solved  for  AY,  it  is  possible  to  use  equation  2.4 
to  calculate  c(y-p).  Therefore,  u(yj)  becomes  a  known  quantity.  An  alternate  method 
of  determining  c(y-j-)  would  be  by  the  use  of  Snell's  law,  or  equation  2.1,  since  P(>'r)» 
P(\'j)  and  c(yj^)  are  all  known. 

Similarly  [Ref.  1:  p.  233], 

kZR  =  kZT  <2-82) 

and,  as  a  result, 


w(yT)  =  T^7w(yR).  (2.83) 

c(yR) 

Referring  to  Figure  2.1  and  utilizing  equation  2.81  it  can  be  seen  that 


c(vn)  HDLTX 

u(yn)  =  cos  a(yn)  =  -^-  u(vR)  =  .  (2.84) 

0  °         c(yR)        K        HRLOS 

Therefore,  since  u(y0)  is  known  from  substituting  y0  for  yy  in  equation  2.81,  the  value 
of  HDLTX  is  given  by 


HDLTX  =  u(y0)HRLOS  .  (2.85) 

Now  that  u(yQ)  and  v(yQ)  are  known  from  equation  2.84  and  equation  2.68,  it  is 
possible  to  find  w(yQ)  by  use  of  the  fact  that  [Ref.  1:  p.  224] 

w2(yn)  =   1  -  u2(yn)  -  v2(yn)  .  (2.86) 
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Again  utilizing  Figure  2.1  and  the  fact  that 
w(y0)  =  cos  y(y0)  ,  (2.87) 

we  see  that 


HDLTZ 

w(yn)  =  .  (2.88) 

0         HRLOS 


Therefore, 

HDLTZ  =  w(yQ)  HRLOS  .  (2.89) 

Figure  2.4  shows  the  geometry  of  Figure  2.1  as  seen  by  looking  down  into  the 
XZ  plane  from  above  the  transmitter's  depth.  From  Figure  2.4  the  relationships 
between  AZ,  AX,  and  AR  may  be  derived. 

The  angle  6  in  Figure  2.4  can  be  found  from 


HDLTX 

tan  6  =  .  (2.90) 

HDLTZ 


Therefore, 

6  =  tan_1(HDLTX/HDLTZ)  .  (2.91) 

Substituting  equations  2.85  and  2.S9  into  equation  2.91  results  in 

6  =  tan-1([u(y0)HRLOS]/[w(y0)HRLOS]}  (2.92) 

so  that 

6  =  tan-'[u(y0Vw(y0)]  .  (2.93) 
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Ax 


Figure  2.4    Topview  of  Geometry. 

For  an  inhomogencous  medium  with  a  sound-speed  profile  that  is  a  function  of 
depth  only,  it  can  be  shown  that  [Ref.  1:  p.  232] 


u(y) 
w(y) 


=  constant 


(2.94) 


Therefore, 


6  =  tan'VypVw^)] 


(2.95) 


where  u(y«)  and  w(yR)  are  available  from  the  frequency  domain  adaptive  bcamformer. 
From  Figure  2.4,  AZ  and  AX  arc  given  by 


AZ  =  AR  cos  6 


(2.90) 


and 
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AX  =  AR  cos  6  .  (2.97) 

Referring  once  again  to  Figure  2.1,  RLOS  is  given  by 

RLOS  =  (AX2  +  AY2  +  AZ2)1/2  (2.98) 

or,  since 

AR2  =  AX2  +  AZ2  ,  (2.99) 


RLOS  =  (AR2  +  AY2)1/2  .  (2.100) 

Finally,  PLOS  can  be  determined  by  using 

PLOS  =  cos_1(AY/RLOS)  .  (2.101) 

The  equations  presented  in  this  section  comprise  the  theory  used  to  develop  the 
model-based  localization  algorithm.  By  the  use  of  ray  acoustics  and  the  assumption 
that  the  model-based  phase  weight  is  known,  a  closed  form  solution  is  possible  for  the 
localization  problem.  Obviously,  the  solution's  accuracy  depends  on  a  ship's  ability  to 
correctly  measure  the  sound-speed  profile  and  the  effects  of  any  other  local  sonar 
conditions,  such  as  shallow  depths  and  the  presence  of  biologies.  However,  in  the  open 
ocean,  when  the  transmitter  and  receiver  are  located  in  the  same  gradient  of  the  sound- 
speed  profile,  a  reasonably  accurate  solution  is  possible.  There  are  some  limitations 
involved  with  the  use  of  ray  acoustics  and  the  model-based  phase  weights.  These 
limitations  will  be  discussed  in  the  next  section. 

C.       LIMITATIONS  OF  RAY  ACOUSTICS  SOLUTION 

1.  Turning  Points 

A  turning  point  is  that  position  along  a  ray  path  propagating  through  an 
inhomogeneous  medium  at  which  the  angle  of  propagation  measured  with  respect  to 
the  positive  Y  axis,  P(y),  is  equal  to  90  degrees.  At  this  point  the  origination  of  the  ray 
path  becomes  ambiguous  to  a  receiver  using  the  localization  technique  described  in  this 
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thesis,  because  there  is  no  way  of  knowing  how  many  turning  points  the  acoustic 
signal  has  passed  through.  The  turning  point  will  cause  a  transmitter  that  is  below 
(above)  the  receiver  to  appear  to  be  above  (below)  the  receiver.  Figure  2.5  illustrates 
these  two  possibilities.  In  the  case  of  receiver  one  in  Figure  2.5,  a  turning  point  has 
occurred  between  the  transmitter's  location  and  the  receiver's  location.  The  theory 
presented  in  this  section  would  result  in  a  calculated  line-of-sight  similar  to  that  shown 
in  Figure  2.5.  The  acoustic  signal  passes  through  two  turning  points  prior  to  reaching 
receiver  two,  and  the  resulting  line-of-sight  calculation  would  indicate  that  the 
transmitter  is  at  a  depth  below  receiver  two. 
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Figure  2.5     Turning  Point  Ambiguity. 

The  turning  point  ambiguity  problem  is  not  necessarily  very  restrictive, 
depending  on  local  sonar  conditions.  Table  1  lists  the  location  of  turning  points  in 
terms  of  AY  and  AR  between  the  transmitter  and  receive  array.  The  values  in  Table  1 
were  calculated  by  assuming  the  values  for  P(y0),  c(y0),  and  g  shown  in  Table  1,  and 
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then  using  equations  2.11  and  2.17  with  p(yj^)  =  90°.  These  results  show  that  for 
most  angles  of  transmission  the  floor  of  the  ocean  would  be  reached  prior  to  the  signal 
reaching  a  turning  point.  Even  at  angles  of  transmission  greater  than  60  degrees  the 
ranges  to  a  turning  point  are  quite  large. 


TABLE  1 

DEPTH  AND  RANGE  TO  TURNING  POINTS  FOR  A  POSITIVE 

GRADIENT 

P(y0) 

AY  (km) 

AR  (km) 

10° 

412.913 

492.043 

20° 

166.935 

238.343 

30° 

86.765 

150.276 

40° 

48.241 

103.424 

50° 

35.921 

86.765 

60° 

13.449 

50.063 

70° 

5.570 

31.582 

80° 

1.336 

15.271 

85° 

0.331 

7.592 

c(y0)  - 

1475 

m/sec 

g  =  0.017 

sec"1 

If  the  transmitter  and  receive  array  are  located  in  the  negative  gradient 
portion  of  the  sound-speed  profile  as  shown  in  Figure  2.5,  the  situation  becomes  much 
more  restrictive.  Here  the  transmitter  must  transmit  in  the  upward  direction  to  reach  a 
turning  point,  as  opposed  to  the  downward  transmission  assumed  in  Table  1.  Table  2 
contains  the  results  of  calculations  for  the  turning  points  in  this  region.  In  this  case, 
the  angles  were  only  varied  from  91  degrees  to  100.8  degrees  in  order  to  place  the 
turning  point  within  the  negative  portion  of  the  sound-speed  profile  of  Figure  2.5. 
Even  with  the  higher  magnitude  gradient  used  in  Table  2,  ranges  of  several  thousand 
meters  are  achievable  prior  to  the  turning  point.  Note  that  all  distances  in  Table  2  are 
in  meters,  whereas  those  listed  in  Table  1  are  in  kilometers. 
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TABLE  2 

DEPTH  AND  RANGE  TO  TURNING  POINTS  FOR  A  NEGATIVE 

GRADIENT 

P(y0) 

AY  (m) 

AR(m) 

91.0° 

-7.601 

870.982 

93.0° 

-68.478 

2615.070 

95.0° 

-190.604 

4365.554 

97.0° 

-374.729 

6126.767 

99.0° 

-621.991 

7903.148 

100.0° 

-769.765 

8798.454 

100.2° 

-801.279 

8978.159 

100.4° 

-833.451 

9158.091 

100.6° 

-866.283 

9338.253 

100.8° 

-899.777 

9518.650 

c(y0)  =  1475 

m 

sec 

g  =  -0.02956 

sec" 

2.  Changes  in  Sound-Speed  Profile 

The  transmitter  and  receiver  must  be  in  the  same  gradient  of  the  sound-speed 
profile  for  the  theory  presented  in  this  thesis  to  work.  If  the  transmitter  and  receiver 
were  located  in  different  gradients  of  the  sound-speed  profile,  a  false  location  would  be 
indicated  due  to  the  change  in  local  angle  of  arrival.  This  situation  is  illustrated  in 
Figure  2.6. 

3.  Validity  of  Model-Based  Phase  Weights 

The  development  of  the  model-based  phase  weights  is  based  in  part  on  the 
assumption  presented  by  Ziomek  [Ref.  1:  p. 253]  that  if 

l[n2(y)-l]'v2(y0)|  <<   1.  (2.102) 

then 

kY(y)  *  kY  +  k02[n2(y)  -  1]  (2kY)  (2.103) 
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Figure  2.6    Changing  Sound-Speed  Gradient. 


n(y)  = 


c(y0) 

c(yR) 


(2.104) 


and 


kY  =  k(>'o)v(>0)  = 


2rtf 

c(y..) 


v(yn)  • 


(2.105) 


For  some  cases,  such  as  P(y0)  approaching  90  degrees,  v(yQ)  becomes  very 
small,  resulting  in  the  criteria  of  equation  2.102  being  violated.  In  these  instances  the 
model-based  phase  weights  can  no  longer  be  considered  valid.  Computations  were 
performed  prior  to  running  the  test  cases  presented  in  this  thesis  to  ensure  that  test 
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cases  which  violate  equation  2.102  were  identified  and  not  misrepresented  as  valid  test 
cases. 

In  addition,  the  WKB  approximation,  which  is  the  basis  for  the  development 
of   the    model-based    phase    weights,    becomes    invalid    as    ky(y)    approaches    zero 
[Ref.  1:  p.  213].   This  is  the  case  at  a  turning  point. 
4.  Depth  Separation  of  Zero  Meters 

If  AY  =  0.0,  meters  the  angle  of  transmission,  P(y0),  and  the  local  angle  of 
arrival,  P(y^),  must  both  be  equal  to  90  degrees  to  permit  the  receive  array  to  receive 
any  signal  without  that  signal  having  to  pass  through  a  turning  point.  The  algorithm 
fails  here  due  to  its  invalidity  at  turning  points  and,  as  can  be  observed  in  equation 
2.17,  because  AR  would  always  be  computed  as  zero.  Obviously,  a  AY  =  0.0  meters 
does  not  necessarily  imply  that  AR  =  0.0  meters,  since  this  condition  is  normally 
known  as  a  collision. 
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III.  COMPUTER  IMPLEMENTATION  OF  LOCALIZATION  THEORY 

A.       PROGRAM  DESCRIPTION 

The  implementation  of  the  theory  described  in  Chapter  II  was  performed  by 
writing  the  FORTRAN  computer  program  LOCATE.  LOCATE  is  designed  to  operate 
as  a  subroutine  in  the  frequency  domain  adaptive  beamforming  algorithm  developed  by 
Ziomek  and  Chan  [Ref.  2].  LOCATE  contains  one  subroutine,  PLOTER,  which 
creates  plots  of  the  function  described  by  equation  2.64.  The  description  of  LOCATE 
that  follows  demonstrates  the  relationship  between  the  equations  of  Chapter  II  and  the 
flow  diagrams,  however,  the  actual  FORTRAN  statements  are  not  presented.  After 
LOCATE  is  explained,  there  is  a  short  discussion  of  PLOTER.  Section  B  discusses  the 
method  by  which  the  algorithm  was  validated.  Section  C  provides  the  actual  results  as 
compared  to  known  geometries,  and  gives  a  comparison  of  double  precision  versus 
single  precision  results. 

1.   Program  LOCATE 

The  program  LOCATE  uses  as  inputs  the  estimated  direction  cosines  for  local 
angles  of  arrival,  model-based  phase  weights,  and  knowledge  of  the  local  sound-speed 
profile  to  determine  AZ,  cross-range  (AX),  depth  separation  (AY),  and  the  line-of-sight 
range  (RLOS)  to  the  transmitter.  Also,  elevation/depression  angle  and  azimuthal  angle 
to  the  transmitter  are  provided  by  LOCATE. 

The  elevation/depression  angle,  as  shown  in  Figure  3.1,  is  defined  as  the 
minimum  angle  between  the  receive  planar  sonar  array's  XZ  plane  and  the  line-of-sight 
between  the  transmitter  and  the  receive  array.  The  elevation/ depression  angle  is 
defined  to  be  positive  (elevation)  if  the  transmitter's  depth  is  less  than  the  receiver's 
depth.  If  the  transmitter  is  at  a  greater  depth  than  the  receive  array  the 
elevation/depression  angle  is  negative  (depression).  Therefore  the  elevation/depression 
angle  ranges  in  value  from  -90  degrees  to  +  90  degrees. 

The  azimuthal  angle,  as  shown  in  Figure  3.2,  is  defined  as  the  minimum  angle 
between  the  receive  planar  sonar  array's  Z  axis  and  the  line-of-sight  between  the 
transmitter  and  receive  array,  in  the  receive  array's  XZ  plane.  The  azimuthal  angle 
then  ranges  from  +  ISO  degrees  to  0  degrees  for  positive  AX  and  from  0  degrees  to 
-180  degrees  for  negative  AX. 

The  inputs  to  the  program  LOCATE  are  defined  as  follows: 
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Figure  3.1     Elevation; Depression  Angle. 

•  UVR,  VYR,  WYR        estimates   of  direction  cosines   u(y^),   v(y^),   and  w(y^), 

respectively,  as  calculated  by  the  frequency  domain 
adaptive  beamformcr. 

model-based  phase  weights. 

carrier  frequency  of  the  received  electrical  signal. 

fundamental  frequency  of  the  finite  Fourier  scries 
representation  of  the  complex  envelope  of  the  received 
electrical  signal. 

gradient  of  local  sound-speed  profile. 

speed  of  sound  at  receive  array  depth  y^. 

total  number  of  receive  elements  along  the  receive  array's  Y 
axis. 

•  QPRIME,  QTOTAL    parameters  used  to  determine  which  harmonic  is  to  be  used 

in  current  calculations. 

•  NPRIME  parameter  used  to  determine  which  element's  phase  weight 

to  use. 

All  the  inputs  are  currently  available  from  the  frequency  domain  adaptive 
bcamforming  algorithm  described  by  Ziomck.  and  Chan  [Ref.  2],  with  the  exception  of 
PHI.    Figures  3.3  and  3.4  illustrates  the  flow  of  the  program  LOCATE. 


•  Fill 

• FREQC 

•  FO 


•  G 

•  CYR 

•  N  TOTAL 
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Figure  3.2    Azimuthal  Angle. 

The  bcaniforming  algorithm  is  written  in  single  precision  FORTRAN. 
However,  the  program  LOCATE  must  operate  in  double  precision  to  enable  it  to 
develop  accurate  roots  for  equation  2.63.  Therefore,  the  values  passed  to  LOCATE 
from  the  adaptive  beamforming  algorithm  must  be  converted  to  double  precision, 
cither  in  LOCATE,  or  before  they  are  sent  to  LOCATE.  In  this  thesis,  all  values 
passed  to  LOCATE  were  double  precision  values.  For  testing  purposes,  only  the 
portions  of  the  adaptive  beamforming  program  which  develop  values  required  by 
LOCATE  were  used,  along  with  a  program  entitled  SOUND  RAY,  which  generates  the 
true  problem  geometry.  The  reasons  for  the  use  of  double  precision  and  the  support 
programs  used  in  testing  LOCATE  are  further  described  in  Section  III.B.l. 

Once  the  program  LOCATE  is  entered,  a  loop  parameter 
QTEMP  =  1,  QTOTAL  is  established.  From  QTEMP,  an  index  Q  for  the  harmonic 
of  interest  is  chosen.  This  value  Q  is  then  used  to  determine  the  frequency  F  that  will 
be  used  for  further  computations. 

To  calculate  the  ray  parameter  SMB  the  local  angle  of  arrival,  P(y^),  is  first 
found  by  the  arc  cosine  of  VYR.    Then  SMB  is  calculated  by  equation  2.2,  using  CYR 
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Fieure  3.3     Program  LOCATE  Flowchart. 
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Figure  3.4     Program  LOCATE  Flowchart. 
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and  P(yj^).  The  value  of  G  is  passed  to  LOCATE  by  the  adaptive  beamformer,  so 
SMA,  the  radius  of  curvature,  is  now  found  by  using  equation  2.12.  At  this  point, 
equations  2.55  through  2.58  are  utilized  to  determine  the  coefficients  A,  B,  C,  and  D. 
These  coefficients  are  in  turn  used  to  find  the  coefficients  of  equation  2.63,  which  are 
stored  in  an  array  called  COEFF. 

To  determine  the  roots  of  equation  2.63,  the  double  precision  IMSL 
subroutine  ZRPOLY  is  called,  using  the  array  COEFF  as  the  input.  ZRPOLY  returns 
complex  roots  for  equation  2.63  in  an  array  called  LAMBDA.  In  all  the  test  cases  that 
were  run,  the  four  roots  in  LAMBDA  always  consisted  of  two  real  roots  and  two 
complex  roots.  As  a  check  of  the  validity  of  the  roots,  the  value  of  F(x,y)  from 
equation  2.64  was  calculated.  A  graph  of  the  function  F(x,y),  such  as  that  shown  in 
Figure  3.5,  was  used  to  determine  whether  to  use  +(1  -  x2)1,2  or  -(1  -  x2)1  2  in  this 
computation  of  F(x,y). 

The  graphs  indicated  that  for  positive  values  of  VYR  the  real  roots  are 
associated  with  the  -+-(1  -  x2)1  2  term,  while  the  complex  roots  are  associated  with  the 
-(1  -x)'  term.  This  can  be  seen  in  Figure  3.5  where  the  curve  associated  with 
-(1  -  x*")1  does  not  cross  the  F(x,y)  =  0  line.  The  graph  in  Figure  3.5  only  shows  a 
small  portion  of  the  X  axis.  Test  runs  demonstrated  that  F(x,y)  increases  as  x  varies 
from  the  x  value  corresponding  to  the  minimum  value  of  F(x,y),  in  both  the  positive 
and  negative  X  directions  over  the  range  0  ^  x  ^  1.  Therefore,  the  graphs  were 
expanded  in  the  region  close  to  the  minimum  of  F(x,y)  to  provide  better  resolution. 

To  continue  with  the  calculations,  one  of  the  four  roots  must  be  selected  as 
the  value  x  of  equation  2.47.  No  logic  in  the  theory  section,  however,  provides  any 
basis  for  a  decision  as  to  which  root  is  correct.  The  complex  roots  were  disregarded 
because  they  cannot  equate  to  x  in  equation  2.47.  In  order  to  determine  a  relationship 
which  would  allow  programming  logic  to  select  the  correct  root  from  the  two  real  roots 
found  bv  ZRPOLY,  numerous  test  cases  with  known  transmitter  and  receive  arrav 
locations  were  run  using  the  four  possible  geometries  allowed  by  the  constraints  listed 
on  page  six  of  this  thesis.   These  geometries  are: 

1.  transmitter  above  receive  array,  0°    ^  P(Yq)  <    90°,   G  >  0. 

2.  transmitter  below  receive  array,  90°  <  P(y0)  <   180°,  G  >  0. 

3.  transmitter  above  receive  array,  0°    <  P(y0)  <    90°,   G  <  0. 

4.  transmitter  below  receive  array,  90°  <  P(y0)  <  180°,  G  <  0. 
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Figures  3.5  through  3.8  correspond  to  each  of  the  four  types  of  geometries 
listed  above.  Analysis  of  these  graphs,  along  with  knowledge  of  the  true  geometry  for 
each  case,  determined  that  if  the  product  G  *  VYR,  designated  RTSLCT  (root 
selection)  in  Figure  3.4,  is  negative,  the  largest  real  root  is  the  correct  value  of  x.  If 
RTSLCT  is  positive,  then  the  smallest  root  is  the  correct  root.  All  test  cases  which 
were  subsequently  run  using  this  root  selection  logic  resulted  in  the  correct  localization 
of  the  target. 

The  root  selected  corresponds  to  x  in  equation  2.47  and  is  next  used  to 
calculate  DELTAY  (AY)  and  DELTAR  (AR),  using  equations  2.11  and  2.17, 
respectively.  From  DELTAY  and  DELTAR  ,  RLOS  is  computed  from  equation  2.100. 
The  azimuthal  angle  is  calculated  next  by  equation  2.95,  since  UYR  and  WYR  are 
known  from  the  adaptive  beamforming  algorithm.  Now  DELTAZ  (AZ)  and  DELTAX 
(AX)  may  be  computed  from  equations  2.96  and  2.97,  respectively. 

The  elevation/depression  angle  is  the  last  value  to  be  computed.  This  is  done 
by  using  equation  2.101,  which  provides  pLOS.  The  angle  PLOS  is  then  converted  to 
the  elevation/depression  angle  by  equation  3.1. 

ELEVDEP  =  90°  -  PLOS  (3.1) 

This  elevation/depression  angle  is  more  useful  than  PLOS  to  personnel 
onboard  ship  because  it  provides  a  target  location  that  is  referenced  to  own  ship's 
horizontal  plane.  Note  that  computing  ELEVDEP  in  this  manner  results  in  a  negative 
value  if  PLOS  >  90°,  which  indicates  that  the  transmitter  is  below  the  receive  array, 
and  a  positive  value  when  pLOS  <90°,  which  implies  that  the  transmitter  is  above  the 
receive  array. 

Program  LOCATE  next  calls  the  subroutine  PLOTER,  if  desired,  to  generate 
a  plot  of  F(x,y)  such  as  that  shown  in  Figure  3.5.  Once  the  graphing  subroutine  is 
completed,  LOCATE  checks  the  index  Q  to  determine  if  the  required  number  of 
harmonics  have  been  evaluated,  and  proceeds  to  process  another  harmonic  if  this  has 
not  been  done.  Otherwise,  the  program  returns  to  the  adaptive  beamforming  program. 
2.  Subprogram  PLOTER 

The  purpose  of  subprogram  PLOTER  is  to  provide  a  graphic  representation 
of  the  roots  which  the  IMSL  subroutine  ZRPOLY  calculates.  The  inputs  to  this 
subprogram  are: 

•  A,  B,  C,  D  coefficients  for  equation  2.64. 
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Figure  3.5     F(x.y)  for  Geometry  1. 
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Figure  3.6     F(x,y)  for  Geometry  2. 
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Figure  3.7     F(x,y)  for  Geometry  3. 
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Figure  3.8     F(x,y)  for  Geometry  4. 
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•  G  gradient  of  the  local  sound-speed  profile. 

•  DELTAX,  DELTAY,  DELTAZ     cross-range,       depth,       and       Z       coordinate 

separations      calculated      by      the      program 
LOCATE. 


pLOS 


the      line-of-sight     angle     as     calculated      by 
LOCATE. 


The  values  of  G,  DELTAX,  DELTAY,  and  DELTAZ  are  printed  out  on  the  graph  as 
G,  AX,  AY,  and  AZ,  respectively,  to  provide  a  means  of  identifying  the  geometry  of 
the  case  corresponding  to  each  graph. 

The  values  A,  B,  C,  D,  and  G  are  converted  to  single  precision  values  prior  to 
being  passed  from  LOCATE  to  PLOTER,  because  PLOTER  was  written  using 
DISSPLA  which  operates  only  in  single  precision.  Due  to  the  single  precision  accuracy 
of  DISSPLA,  plots  made  by  PLOTER  are  not  accurate  enough  to  determine  the  roots 
of  equation  2.64.  However  the  plots  do  show  approximately  where  the  roots  occur. 
Figure  3.9  illustrates  the  flow  of  the  subroutine  PLOTER. 


START 


DETERMINE  X  COORDINATE  FOR 
MINIMUM  F(X.Y)  =  XMIN 


I 


COMPUTE  F(X,Y)  FOR 
XMIN  -  0.025  <  X  <  XMIN  +  0.025 


I 


PLOT  F(X,Y)  VERSUS  X 


I 


RETURN 


Figure  3.9     Subprogram  PLOTER  Flowchart. 
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The  subroutine  PLOTER  first  computes  the  minimum  value  of  F(x,y)  in  the 
interval  0  <  x  <  1  by  incrementing  x  by  0.1  units.  This  minimum,  XMIN,  is  then 
used  as  the  center  of  the  plot,  with  XMIN  -  0.025  and  XMIN  +  0.025  as  the  lower 
and  upper  bounds  of  the  graph.  If  XMIN  +  0.025  >  1.0  the  plot  is  centered  at  0.975 
to  avoid  having  the  computer  attempt  to  calculate  the  square  root  of  a  negative  value 
of  1  -  x2  in  equation  2.64.  After  the  plot  is  completed,  PLOTER  returns  to  the 
program  LOCATE. 

B.       ALGORITHM  VALIDATION 

1 .  Generation  of  Received  Signals 

The  inputs  listed  in  Section  A  of  this  chapter  for  the  program  LOCATE  were 
generated  through  the  use  of  two  programs.  The  first  program  is  titled  SOUNDRAY 
and  was  written  by  Professor  L.  J.  Ziomek  at  the  U.  S.  Naval  Postgraduate  School, 
Monterey,  California,  in  1987.  The  second  program  is  the  subroutine  PHSWGT 
developed  by  Ziomek  and  Blount  [Ref.  7].  SOUNDRAY  utilizes  ray  acoustics  and 
geometry  to  develop  feasible  geometries  for  calculations  of  local  angles  of  arrival  of 
acoustic  signals.  The  inputs  to  SOUNDRAY  are  the  X,  Y,  and  Z  coordinates  of  the 
transmitter,  the  X  and  Y  coordinates  of  the  receive  array,  the  initial  angle  of 
propagation,  P(y-p),  and  information  describing  the  local  sound-speed  profile. 
SOUNDRAY  then  uses  equation  2.1  to  determine  P(y^)  and  equation  2.15  to  calculate 
AR.  From  this  point,  geometry  alone  allows  calculation  of  the  RLOS  and  PLOS,  from 
equations  2.98  and  2.101,  respectively,  and 

AZ  =  (RLOS2  -  AX2  -  AY2)1*'2  .  (3.2) 

In  addition,  SOUNDRAY  calculates  the  inputs  for  the  subroutine  PHSWGT 
and  the  estimates  (in  this  case  exact  values)  of  direction  cosines  for  the  acoustic  signal 
arriving  at  the  receive  array.  SOUNDRAY  determines  the  exact  problem  geometry, 
independent  of  the  model-based  phase  weights,  thereby  providing  the  standard  by 
which  to  judge  the  solutions  generated  by  the  program  LOCATE. 

2.  Test  Case  Results 

a.   Double  Precision  LOCA  TE  versus  True  Geometry 

As  stated  previously,  there  are  four  basic  geometries  that  the  program 
LOCATE  is  designed  to  handle.   These  four  geometries  may  be  summarized  as: 
1.      +AY,  0°  <  P(yQ)  <    90°,  G  >  0. 
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2.  -AY,  90°  <  P(yQ)  <  ISO0,  G  >  0. 

3.  +  AY,  0°  ^  P(yQ)  <    90°,  G  <  0. 

4.  -AY,  90°  <  P(yQ)  <   180°,  G  <  0. 

Other  variations  on  these  geometries  are  possible  by  using  -AX  and  -AZ, 
but,  because  the  sound-speed  profile  is  assumed  to  be  a  function  of  depth  only,  the 
plane-wave  field  will  propagate  in  a  plane  which  is  normal  to  the  XZ  plane  [Ref.  1:  p. 
234].  The  result  is  that  variations  using  -AX  and  -AZ  merely  change  the  sign  of  the 
solutions  and  not  the  magnitude.  LOCATE  was  written  to  accommodate  -AX  and 
-AZ.  However,  for  this  discussion,  it  is  sufficient  to  deal  with  +AX  and  +AZ  and 
realize  that  only  the  sign  of  the  answer  is  different  when  negative  quantities  are  used. 

Tables  3,  4,  5,  and  6  represent  results  from  the  four  geometries  mentioned 
above.  The  sound-speed  profile  of  Figure  2.2  was  used  in  these  computations.  The 
value  of  AY  for  each  table  was  maintained  constant  and  this  necessitated  the  altering 
of  AX  depending  on  the  angle  P(y0)  used.  If  P(y0)  was  close  to  0  degrees  or 
180  degrees,  a  smaller  AX  was  required  than  for  angles  near  90  degrees.  This  is  due  to 
the  fact  that  at  angles  near  0  degrees  or  180  degrees,  the  plane-wave  field  reaches  depth 
y^  in  a  much  shorter  AR  than  when  P(y0)  is  near  90  degrees.  Since  from  equation 
2.99 

AR2  =  AX2  +  AZ2  ,  (2.99) 

AX  had  to  be  kept  sufficiently  small  to  maintain  AZ  >  0,  because  we  are  working  with 
cases  of  positive  AX  and  AZ. 

As  can  be  seen  in  Tables  3  through  6,  the  program  LOCATE  provides 
excellent  results.  The  slight  errors  that  are  present  are  due  mainly  to  roundoff  error 
occurring  in  the  root  finding  subroutine  ZRPOLY.  Note  that  the  constraints 
concerning  turning  points  have  all  been  observed  in  these  results.  The  maximum  error 
for  any  range  calculated  by  LOCATE  in  these  cases  was  0.345  meters.  The  angles 
calculated  by  LOCATE  are  not  presented  in  tabular  form  because  they  were  all 
accurate  to  four  significant  digits  when  compared  to  the  true  solutions. 

Some  of  the  results  in  Tables  3  through  6  appear  to  be  exact.  This  is  not 
actually  the  case  because  the  values  in  these  tables  were  all  rounded  to  the  third 
decimal  place.  In  no  instance  were  the  results  of  LOCATE  exactly  equal  to  the  true 
solution,  however,  in  many  instances,  the  difference  was  in  the  fourth  or  fifth  decimal 
place. 

51 


TABLE  3 

LOCATE  VERSUS  TRUE  GEOMETRY  (GEOMETRY  1) 

P(y0) 

AX(m) 

AY(m) 

AZ( 

m) 

T 

L 

T 

L 

T 

L 

10° 

50.000 

50.037 

300.000 

300.000 

17.434 

17.449 

15° 

50.000 

50.017 

300.000 

299.999 

63.096 

63.118 

20° 

100.000 

100.033 

300.000 

299.999 

44.287 

44.303 

25° 

100.000 

100.017 

300.000 

299.999 

98.212 

98.229 

30° 

100.000 

100.015 

300.000 

300.000 

141.879 

141.900 

35° 

100.000 

100.017 

300.000 

300.000 

185.308 

185.340 

40° 

100.000 

100.010 

300.000 

299.999 

231.794 

231.819 

45° 

300.000 

300.037 

300.000 

300.000 

24.554 

24.559 

50° 

300.000 

300.020 

300.000 

300.000 

197.192 

197.206 

55° 

300.000 

300.029 

300.000 

300.000 

308.993 

309.023 

60° 

500.000 

500.131 

300.000 

300.001 

153.761 

153.801 

65° 

500.000 

500.121 

300.000 

300.001 

414.580 

414.680 

70° 

500.000 

500.049 

300.000 

299.976 

670.746 

670.811 

75° 

500.000 

500.064 

300.000 

300.001 

1035.436 

1035.569 

80° 

500.000 

500.038 

300.000 

300.001 

1741.051 

1741.185 

85° 

500.000 

500.019 

300.000 

300.000 

5226.883 

5227.089 

T  =  true  solution 

L  = 

LOCATE  calculation 

G  =   +0.017  sec"1 

b.   Errors  as  a  Function  of  Angle  of  Transmission  and j or  Depth  Separation 

(1)  Depth  Separation.  Figure  3.10  shows  the  error  in  RLOS  as  the  depth 
separation  between  the  transmitter  and  the  receive  array  increases,  with  P(y0)  constant. 
There  does  not  seem  to  be  any  relation  between  the  error  and  the  depth  separation. 
The  error  appears  to  be  mainly  caused  by  roundoff. 

(2)  Transmission  Angle  and',  or  Depth  Separation.  Figure  3.11  shows  the 
error  in  RLOS  as  the  angle  of  transmission  changes  for  four  different  depth 
separations.  Again,  it  is  readily  observed  that  the  depth  separation  has  little  effect  on 
the   size  of  the  error. 
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TABLE  4 

LOCATE  VERSUS  TRUE  GEOMETRY 

(GEOMETRY  2) 

P(y0) 

AX(m) 

AY(m) 

AZ( 

m) 

T 

L 

T 

L 

T 

L 

95° 

500.000 

500.060 

-300.000 

-299.993 

2907.664 

2908.009 

100° 

500.000 

500.056 

-300.000 

-299.985 

1536.780 

1536.952 

105° 

500.000 

500.013 

-300.000 

-299.976 

971.804 

971.831 

110° 

500.000 

500.037 

-300.000 

-300.000 

640.734 

640.782 

115° 

500.000 

500.052 

-300.000 

-300.000 

395.290 

395.332 

120° 

500.000 

500.062 

-300.000 

-300.000 

128.008 

128.024 

125° 

400.000 

400.027 

-300.000 

-300.000 

147.308 

147.318 

130° 

300.000 

300.129 

-300.000 

-300.000 

191.553 

191.637 

135° 

200.000 

200.065 

-300.000 

-300.000 

222.138 

111.211 

140° 

200.000 

200.022 

-300.000 

-300.000 

151.643 

151.660 

145° 

200.000 

200.051 

-300.000 

-300.000 

62.335 

62.353 

150° 

100.000 

100.072 

-300.000 

-300.000 

140.799 

140.901 

155° 

100.000 

100.062 

-300.000 

-300.000 

97.297 

97.358 

160° 

100.000 

100.066 

-300.000 

-300.000 

43.151 

43.182 

165° 

50.000 

50.047 

-300.000 

-300.000 

62.552 

62.722 

170° 

50.000 

50.063 

-300.000 

-300.000 

16.779 

16.804 

T  =  true  solution 

L  = 

LOCATE  calculation 

G  =   +0.017  sec"1 

The  error  does  increase  as  the  angle  P(y0)  is  increased  above  about 
60  degrees.  This  increase  can  be  attributed  to  the  behavior  of  the  sine  and  cosine 
functions.  Figure  3.12  shows  how  the  sine  and  cosine  functions  behave  between  0  and 
90  degrees.  Above  about  60  degrees,  the  slope  of  the  sine  function  is  less  than 
0.01  degrees"1  so  that  small  changes  in  the  sine  cause  large  differences  in  the  angle  P(y). 
Also,  in  this  region  the  magnitude  of  the  slope  of  the  cosine  function  is  near  its 
maximum.   Small  changes  in  the  angle  P(y)  create  large  differences  in  the  cosine. 
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TABLE  5 

LOCATE  VERSUS  TRUE  GEOMETRY  (GEOMETRY  3) 

P(y0) 

AX(m) 

AY(m) 

AZ< 

m) 

T 

L 

T 

L 

T 

L 

10° 

50.000 

49.995 

500.000 

500.000 

72.077 

72.071 

15° 

100.000 

99.994 

500.000 

500.000 

88.099 

88.094 

20° 

100.000 

99.995 

500.000 

500.000 

150.837 

150.830 

25° 

100.000 

100.000 

500.000 

500.000 

209.070 

209.067 

30° 

100.000 

99.997 

500.000 

500.000 

268.783 

268.777 

35° 

300.000 

300.004 

500.000 

500.000 

175.431 

175.434 

40° 

300.000 

300.001 

500.000 

500.000 

288.241 

288.243 

45° 

300.000 

300.009 

500.000 

500.000 

393.838 

393.850 

50° 

500.000 

500.008 

500.000 

499.999 

310.993 

310.999 

55° 

500.000 

500.006 

500.000 

499.999 

494.933 

494.940 

60° 

500.000 

499.979 

500.000 

500.000 

686.650 

686.620 

65° 

500.000 

499.976 

500.000 

500.001 

916.323 

916.279 

70° 

500.000 

499.997 

500.000 

500.000 

1221.188 

1221.181 

75° 

500.000 

499.995 

500.000 

500.000 

1671.182 

1671.169 

80° 

500.000 

500.002 

500.000 

499.999 

2426.104 

2426.112 

85° 

500.000 

500.004 

500.000 

499.998 

3902.854 

3902.891 

T  =  true  solution 

L  = 

LOCATE  calculation 

G  =  -0.02956  sec" 

1 

To  find  AY,  equation  2.11  uses  the  roots  of  equation  2.63  as 
determined  by  ZRPOLY.  These  roots  correspond  to  sin  P(y0).  The  root  contains 
some  small  errors  due  to  roundoff  which  is  borne  out  by  the  fact  that  the  values  of  AY 
in  Tables  3  through  6  contain  errors  on  the  order  of  10"3  meters.  To  find  AR  by  using 
equation  2.17,  the  arc  sine  of  the  root  must  first  be  calculated.  This  amplifies  any  error 
in  the  root,  especially  when  the  angle  is  greater  then  60  degrees  as  discussed  previously. 
Next,  the  cosine  of  the  arc  sine  of  the  root  is  computed,  which  further  amplifies  the 
error. 
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TABLE  6 

LOCATE  VERSUS  TRUE  GEOMETRY  (GEOMETRY  4) 

P(y0) 

AX(m) 

AY  (m) 

AZ( 

m) 

T 

L 

T 

L 

T 

L 

100° 

500.000 

499.998 

-500.000 

-500.000 

3539.164 

3539.153 

105° 

500.000 

500.000 

-500.000 

-500.000 

1966.455 

1966.417 

110° 

500.000 

499.983 

-500.000 

-500.000 

1347.650 

1347.606 

115° 

500.000 

500.000 

-500.000 

-500.000 

983.983 

983.992 

120° 

500.000 

499.986 

-500.000 

-500.000 

728.904 

728.884 

125° 

500.000 

499.961 

-500.000 

-500.000 

525.292 

525.252 

130° 

500.000 

499.965 

-500.000 

-500.000 

337.478 

337.455 

135° 

500.000 

499.980 

-500.000 

-500.000 

71.388 

71.389 

140° 

300.000 

299.990 

-500.000 

-500.000 

298.440 

298.432 

145° 

300.000 

299.979 

-500.000 

-500.000 

185.560 

185.548 

150° 

100.000 

99.991 

-500.000 

-500.000 

272.886 

272.863 

155° 

100.000 

99.987 

-500.000 

-500.000 

212.229 

212.201 

160° 

100.000 

99.990 

-500.000 

-500.000 

153.303 

153.288 

165° 

100.000 

99.974 

-500.000 

-500.000 

90.286 

90.264 

170° 

50.000 

49.978 

-500.000 

-500.000 

73.210 

73.181 

T  =  true  solution 

L  = 

LOCATE  calculation 

G  =  -0.02956  sec" 

L 

Therefore,  above  about  60  degrees,  we  see  these  increased  errors 
manifest  themselves  in  the  AR,  AX,  AZ,  and  RLOS  calculations.  Still,  the  errors  seen 
in  Figure  3.11  and  in  Tables  3  through  6  are  insignificant  when  compared  with  the 
ranges  in  question.  The  angles  are  still  accurate  to  four  significant  digits,  and 
consequently,  the  range  errors  remain  small. 

c.    Double  Precision  Versus  Single  Precision  Results 

It  was  found  that  the  single  precision  version  of  ZRPOLY  was  not  accurate 
enough  to  calculate  the  correct  answers.  The  reason  for  this  can  be  seen  in  Table  7 
which  contains  some  single  precision  results  for  comparison  to  double  precision  results. 
ZRPOLY  calculates  the  roots  shown  in  the  two  right  hand  columns  of  Table  7.    Even 
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Figure  3.10     Error  in  RLOS  as  a  Function  of  Depth  Separation. 
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Figure  3.1 1     Error  in  RLOS  as  a  function  of  Transmission  Angle  and' or  Depth  Separation. 
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Figure  3.12     Sine  and  Cosine  for  0  to  90  Degrees. 
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TABLE  7 

DOUBLE  PRECISION  VERSUS  SINGLE  PRECISION  RESULTS 

DP 

SP 

DP 

SP 

DP 

SP 

too) 

P(y0) 

RLOS  (m) 

RLOS  (m) 

Root 

Root 

60.17° 

60.41° 

603.147 

9.264 

0.8660 

0.8689 

65.21° 

65.44° 

715.599 

48.291 

0.9063 

0.9092 

70.28° 

70.58° 

888.832 

28.699 

0.9397 

0.9428 

75.38° 

75.75° 

1188.473 

30.835 

0.9659 

0.9693 

80.60° 

81.18° 

1836.237 

57.808 

0.9912 

0.9881 

86.73° 

88.34° 

5259.514 

350.602 

0.9961 

0.9997 

G  =   +0.017 

sec 

though  the  roots  appear  accurate  to  the  second  significant  digit  in  the  single  precision 
results,  when  dealing  with  sines  and  cosines,  an  error  in  the  third  significant  digit  can 
create  a  fairly  large  error  in  calculating  the  angle  P(y0).  Also,  these  roots  are  multiplied 
by  the  radius  of  curvature,  a,  in  equation  2.11.  This  radius  of  curvature  is  on  the  order 
of  10  meters,  so  small  errors  in  the  roots  will  create  large  errors  in  the  ranges 
calculated.  The  single  precision  results  in  Table  7  are  so  poor  that  they  seem  to  have 
no  relation  to  the  actual  answer.  The  double  precision  results  for  RLOS  in  Table  7  are 
accurate  to  within  0.1  meters  of  the  true  solution. 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  goal  of  this  thesis  was  to  determine  if  an  underwater  acoustic  transmitter  can 
be  localized  using  ray  acoustics,  model-based  phase  weights,  estimates  of  the  local 
angles  of  arrival,  and  knowledge  of  the  local  sound-speed  profile.  As  demonstrated  in 
Chapter  III,  this  goal  is  achievable  and  to  a  high  degree  of  accuracy  depending  on  the 
accuracy  of  the  inputs  to  LOCATE.  There  are  restrictions  on  the  use  of  this 
procedure.  It  appears  that  the  restrictions  do  not  impose  severe  limitations  on  the  use 
of  the  algorithm,  and  in  some  cases  it  may  be  possible  to  overcome  them  altogether. 

All  the  restrictions  basiclly  result  in  a  limitation  on  the  effective  range  of  the 
algorithm.  Even  though  acoustic  signals  may  not  reach  their  initial  turning  points  for 
theoretical  ranges  in  the  tens  or  even  hundreds  of  kilometers,  the  ocean  is  only  about 
11.5  kilometers  deep  at  its  greatest  depth.  Therefore,  the  ranges  shown  in  Table  1  are 
not  realizable  in  some  cases  because  the  signal  will  reach  the  ocean  floor  in  less  range 
than  it  would  take  to  reach  the  turning  point.  Additionally,  underwater  acoustic 
transmitters  are  usually  limited  in  the  depth  to  which  they  may  be  deployed,  so  that  the 
angles  of  transmission  that  are  associated  with  the  greatest  ranges  will  pass  well  below 
the  receive  array  at  any  significant  range.  Still,  the  algorithm  appears  to  be  quite 
useable  in  ranges  of  less  than  10  kilometers.  This  would  be  of  a  great  advantage  in  the 
case  of  a  transmitter  whose  signal  is  of  low  power,  resulting  in  a  short  detection  range. 
In  fact,  the  need  for  an  algorithm  of  this  sort  is  most  critical  when  the  transmitter  is  at 
short  range  and  its  exact  location  and  direction  of  motion  must  be  resolved  rapidly. 

In  some  instances,  the  limitations  due  to  turning  points  may  not  be  of  much 
concern.  For  example,  the  algorithm  might  be  used  for  an  array  located  on  the  ocean 
floor.  In  this  case,  much  longer  ranges  would  be  achievable,  provided  that  the 
transmitter  is  in  the  same  portion  of  the  sound-speed  profile  as  the  receive  array.  The 
algorithm  might  also  be  of  use  in  active  sonar  systems  to  provide  more  accurate  range 
and  depth  information  than  is  currently  available. 

Implementation  of  the  algorithm  must  include  a  very  accurate  root  finding 
technique  as  has  been  discussed.  Due  to  the  sensitivity  of  the  problem  in  regard  to  the 
sine  and  cosine  functions,  the  roots  need  to  be  accurate  to  at  least  three  significant 
figures.  It  was  found  that  this  is  only  possible  through  use  of  a  double  precision  root 
finding  subroutine.    This,  of  course,  causes  the    program  to  run  more  slowly  but, 
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because  the  remainder  of  the  program  can  still  be  written  in  single  precision,  it  is  not  a 
great  hinderance. 

In  the  future  some  areas  requiring  more  study  are: 

•  Develop  a  method  for  obtaining  the  model-based  phase  weights  from  the 
received  signals.  At  present,  phase  weights  are  computed  based  on  received 
signals,  however,  the  phase  weights  in  the  Y  direction  need  to  be  separated  into 
traditional  phase  weights  and  model-based  phase  weights. 

•  Determine  a  method  to  account  for  the  acoustic  signal  passing  through  a 
turning  point  prior  to  reaching  the  receive  array.  This  would  greatly  extend  the 
range  capability  of  the  algorithm. 

•  Develop  methods  to  identify  signals  that  are  transmitted  from  portions  of  the 
sound-speed  profile  other  than  the  gradient  in  which  the  receive  array  is  located. 

•  Investigate  the  practical  applications  of  the  algorithm  in  varying  acoustic 
conditions,  particularly  in  regions  such  as  near  the  Gulf  Stream  where  the 
sound-speed  profile  is  a  function  of  depth  and  range. 
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